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Abstract 

Instabilities in finite difference codes due to the singularity of spherical coor- 
dinates at the center are studied. In typical Numerical Relativity applications, 
standard regularization techniques by themselves do not ensure long term sta- 
bility. A proposal to remedy that problem is presented, which takes advantage 
of redundant quantities introduced in recent hyperbolic formulations of Ein- 
stein's evolution equations. The results are discussed through the example 
case of a boson star, where a significant improvement in the implementation 
of boundary conditions is also presented. 

PACS numbers: 04.25.Dm 
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I. INTRODUCTION 



Spherically symmetric systems have deserved a lot of interest since the early stages of 
Numerical Relativity. In the pioneering work of May and White the evolution of a self- 
gravitating spherical star was modelled by coupling the hydrodynamical equations with two 
first integrals of Einstein field equations: the energy and momentum constraints. These are 
two ordinary differential equations, which in spherical symmetry are enough to determine 
the evolution of the whole metric. The validity of this approach, however, is restricted to 
the spherical case, because in more general situations these first integrals are not enough to 
determine the solution of Einstein field equations. Furthermore, even in the spherical case, 
the treatment of boundary conditions with this system is known to be rather difficult. The 
latter is especially dangerous in systems with outgoing matter or radiation, which need to 
be carefully dealt with. 

The alternative approach is to couple the hydrodynamical equations with the full set of 
Einstein partial differential equations. In that way one must evolve both the hydrodynamical 
quantities and the spacetime metric on the same footing. The resulting formalism can then 
be extended in principle to axially symmetric (2D) or even completely general situations. 
In the spherically symmetric case, however, the major drawback of this approach is not just 
the higher number of metric quantities to compute. It is rather the fact that the spherical 
coordinate system is singular at the origin, even when the spacetime itself is regular. The 
region close to the center of symmetry poses serious problems to numerical codes, which are 
to be overcome by using special techniques like artificial viscosity or algorithm matching. In 
most cases, this is usually done by trial and error and, together with the treatment of the 
outer boundary, becomes the main limiting factor when constructing finite difference codes. 

The purpose of this paper is to present a new way of writing the full set of Einstein 
equations in spherical symmetry (with a matter content), which takes advantage of some 
redundant quantities that appear in recent hyperbolic formulations of Numerical Relativ- 
ity in order to solve this problem. In that way, one will be able to treat the center 
region with exactly the same numerical algorithms as any other region of the numerical grid, 
with a significative reduction of the time devoted to code development. The same formalism 
also allows us to impose more accurate external boundary conditions, related to the physics 
of the problem. The combination of these two improvements (center and boundary treat- 
ment) helps in increasing the algorithm stability and maintaining the required accuracy. As 
far as tested, the resultant code does not show any stability problem, even after very long 
runs. We will work out the description of the method through the case of a complex massive 
scalar field describing a radiating boson star, where we easily reproduce previous results 
from other authors [Q. 

II. FORMALISM 

The object of Numerical Relativity is to solve the Einstein equations 

G^u = SirGTfj^u, (1) 

in an approximate way. This is necessary when we deal with some of the most interesting 
systems in astrophysics: black holes, neutron stars, close binaries, supernovae, etc. 
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G^u is the Einstein tensor associated to a certain spacetime (represented by a metric (7^,^), 
G is the gravity constant and T^j^ is the momentum-energy tensor describing matter in this 
spacetime. 

The metric element can be written, in normal coordinates, as 

ds^ = g^^dx^dx'^ = —a^dt^ + dx'^ dx^ (2) 

where = gij is the space metric and a is the lapse, a gauge degree of freedom. 
The Einstein system (|Tp is constituted by 10 coupled second order partial differential equa- 
tions, when written in terms of the 4-dimensional metric g^v Six of them are evolution 
equations and the four remaining are constraint equations 

We deal with spherical symmetry problems, so the system reduces to a one-dimensional (ID) 
problem. In spherical coordinates, (|^) is written 

— a^dt^ + grrdr"^ + gee dVL (3) 

where dVt = dO"^ + sin? 9 dcj)^. 

The Einstein system can be written in first order form by introducing new variables. For 

simplicity, we will use the vacuum equations to illustrate our procedures. Following we 
will express it as it follows: 

dtgee = gee{Cq%) , (4) 

dtgrr = grriCq\) , (5) 



dtC = C{l/2C[if-l)q\ + 2fq\]) , (6) 

dtLr - dr {C f[q\ + 1/2 q\]) = , (7) 

dtDre'-dr{Cq\)=0, (8) 

dttrD - dr {C{q\ + 2g%)) = , (9) 

dtVr. = 2C [q\ Lr - l/2(g% - g^,) A/] , (10) 

dtiq^r) - dr (2C[3/4 Vr - A./ + U]) = 

-c(l/2 q%{q'r-q\)-^— 
V gee 

+ Lr Dre' + l/^[{q\f - [Dre'f]) , (11) 
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dt{q\) - Or {C[Dj - 1/2 Vr]) = 

-c(l/2 q\iq\-q\) + ^ 
V gee 

- U Dj + l/A[{q\f - [Djf]) 



(12) 



where we have defined = drln{a), Dr/ = drln{grr), Dre^ = drln^geg), trD = Dr/+2Dre , 
C = — = (speed of hght), q^r = ^ind q^g = ^^l^- Also, allowing for the momentum 



constraint, Vr = 2Dre ■ 

Note that we have written the Einstein system as a first order balance law system. The vari- 
able Vr is redundant: it was introduced with the purpose of making the system hyperbolic, 
which allows us to directly deal with diagonalized variables with a well-defined speed and 
sense of propagation (very useful when imposing boundary conditions, as we will see later). 

The function /(a) determines the gauge. We have chosen the gauge degree of freedom 
a to follow the harmonic slicing which is the simplest choice (/ = 1) and good enough 
for singularity-free spacetimes. 



Due to spherical symmetry ggg has a strong geometrical factor: it is proportional to 
when r ^ 0. It follows that both D^e^ and Vr are proportional to 1/r , a singular behaviour 
at the origin. 

Simulations dealing directly with these equations, or equivalent ones, are extremely unstable. 
This is well-known and some alternatives are usually taken. In previous works dealing 
with black holes, the problem was avoided by putting internal boundary conditions at the 
apparent horizon. In this way the grid was cut at the apparent horizon, far enough from 
both the physical and coordinative singularity at r = 0. The problem is completely different 
when we deal with singularity-free systems (stars). There, there is no way of cutting the 
solution before r = without loosing physically relevant information. The coordinative 
singularity at r = turns up as a big danger for accuracy and stability near the origin when 
using the standard system (^[121). 

The first step for dealing with this is to analytically extract the geometrical factors, as far 
as possible, from the equations, so they only deal with the regular part. This is done by 
writing the line element as 



III. GEOMETRICAL FACTORS AND NUMERICAL ERROR 



— a^dt^ + grr dr"^ + r'^gee dVt 



(13) 



so that the system, in our particular case, is rewritten with gso = r'^gee, Dre 
trD = Dr/ + 2Dre^ ■ The new equations read 



e 



dtgee =geeCqg 



(14) 




(15) 
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dttrD - dr (C{q\ + 2q\)) = (16) 

This is yet an improvement because the evolved quantities -0,0^ and gge are now regular 
at the center. But further improvement should be achieved, provided we want to deal with 
a stable code, because rewriting (^-p!^) with (|1^) introduces explicit 1/r factors in the terms 
containing space derivatives (flux terms) and factors show up in the right hand side 
(source terms) of equations (|TTj) and (|T^). We have somehow explicitly segregated the geo- 
metrical term, but it is still not out of the system. Now we can appreciate the kind of error 
that is introduced: both fluxes and sources are singular due to 1/r and 1/r^ terms and exact 
cross-cancellation is needed between them to keep the system regular. However, this is not 
achieved by the finite difference algorithm due to the truncation error when calculating the 
space derivatives of the fluxes. The truncation error when computing a given variable with a 
second order code is proportional to (Ar)^, where Ar is the distance between neighbouring 
points in the grid. The terms 1/r in the fluxes give errors like (Ar)^/r^. Where r — i> 
(r ~ Ar) the code does not converge. 

As a particular instance we tried to solve the boson star ID problem |^ by using these 
modified quantities. The numerical truncation error caused an instability that prevented us 
from obtaining any meaningful result before the code crash. 

These are well-known regularization techniques to handle this problem. In the next section, 
however, we present a new approach which regularizes the equations in a deeper way. 



IV. THE ROLE OF THE MOMENTUM CONSTRAINT 

Our proposal goes beyond the standard regularization procedures and succeeds in remov- 
ing the most dangerous remaining 1/r factors. We take advantage from the way in which the 
momentum constraint was used to build the system (§-^). Remember that the variable K- 
is redundant and it was introduced with the objective of making the system hyperbolic 0. 
As said, in (^[12D the momentum constraint is satisfied if and only iiVr = 2Dre^ ■ We can 
monitorize the level of violation of the momentum constraint by testing this algebraic iden- 
tity through the evolution. 

From the considerations of the preceding section, it will be natural to evolve instead of Vr 
the regular quantity K = '^Dro^. But it turns out to be more convenient to define instead 



Vr^2Dj + - (17) 




which is also regular (provided that Qrr/gee — 1 ^ as r^, which demands smooth initial 
data near the origin) and removes both the 1/r terms in the fluxes and the 1/r^ terms in 
the sources. The singularities have now disappeared from all the evolved variables as well as 
the numerical error caused by the incomplete cross-cancellation between fluxes and sources. 
The modified equations are 



dtVr = C 2q%Lr + {q\ - q%)^ (18) 
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^rd^{— DrO^) + 2LrDre^ 

-hrVr + h^Dj - ^ - trD) ] (19) 
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^ ^ + - {2Lr + AD re' -trD)) (20) 
lb r / 



There is still a point to be taken into account. In eqs. ( p!9| - pO| ) the sources contain terms like 
1/r times trD, Lr or other variables that are radial derivatives of the metric coefficients. 
These terms do not pose any problem as r — * 0, as radial derivatives of any differentiable 
function vanish at the origin. However, in the finite difference approach, we cannot put a 
grid point at r = due to the explicit 1/r factors. We will therefore place the ffist grid 
point at r = Ar/2. With this new way for writing the Einstein equations in spherical 
symmetry we are able to evolve the whole grid, without needing to apply special techniques 
or different algorithms to the center region. The equations are now intrinsically stable and 
can be directly coded. 



V. EXAMPLE: THE BOSON STAR 



We have applied this idea ( p!7D to a Boson Star model g]. Mathematically, the system 
is constituted by the Einstein equations plus the Klein-Gordon ones. It has stable configu- 
rations known as Boson Stars, which determine a mass-radius curve, as in the neutron star 
case, with a local maximum for the star masses. 

The stars located at the left (S-Branch) from this maximum (small radius) are known to be 
stable against small perturbations. The stars at the right (U-Branch) are, on the other hand, 
unstable. An U-Branch star, when perturbed can migrate towards a stable configuration 
in the S-Branch. The migration produces strong periodical oscillations in the star radius 
around the final position, as well as scalar field radiation. 

We tested our equations against this phenomenon. We could not evolve the system with 
■ jT^ ) far enough to measure any relevant physical quantity. With p^p!6| , p!8| - pOD , we 



have been able to reproduce the migration and the complex behaviour of its stationary 
states, with run time just limited by the required accuracy and not by algorithm stability, 
which is not easily achieved with the alternative formulation. We have to point out that 
part of the success comes from the hyperbolic balance law formulation itself, by allowing 
more physically consistent external boundary conditions. In the next section we go further 
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in that point. In figure 1 we show the migration from an unstable configuration to a stable 
one 0. With the original form of (|^- |T2|) we were unable to see nor half an oscillation. 



VI. THE BOUNDARY CONDITIONS 



We illustrate the method for boundary conditions by describing the treatment of the 
Klein-Gordon equations in our example system. 



Let's take a look at the equation for the real part, 
m^0, in tensor form): 



of the Klein- Gordon field {g 



dtpi - driC^i) = 

dtiJl - driCpi) = C 



(^Dr9^ +'^^ Pi- (fei^l - m^grr<Pl 



(21) 
(22) 



where p = (9^0, = -^dtcj). These equations, like the rest of the system, are written in 
balance law form. The flux terms are those that would appear in a fiat space wave equation 
(transport terms). The term in the right hand side is an oscillatory term which comes 
from the m^(f) in the Klein-Gordon equation. The rest of the source terms are due to spherical 
coordinates and the coupling between the Klein-Gordon field and the spacetime metric. 
To put a condition over the transport terms would be quite easy, but the combination with 
the source terms makes the problem more difficult. What one would like to do is to treat 
separately each term of the equation. Fortunately, this is possible in a numerical code, by 
using the 'operator splitting' approach. With the Strang splitting (second order operator 
splitting) each step in the numerical evolution can be decomposed in the following way: 



e{At) = S{^)n{At)S{^) 



(23) 



The whole step, e(At), is split into three sequential processes. The first one is S{f). There 
we evolve during half a step only the 'source' terms in the equations. If we restrict to the 
Klein-Gordon part, we have: 



dtpi 





c 



Dre" + - ) Pi - q%4^i - m grr(j)i 
r , 



(24) 
(25) 



This is an ordinary differential equation system for each point in the net. In this part there 
is no connection between points. Each point evolves by itself. So we do not have to state 
any boundary condition in this substep. 

In the second substep, 72.(Ai:), we only evolve the transport part (during a whole time step), 
so we have 



9iPi-9,(C^i)=0 
- 9,(Cpi) = 



(26) 
(27) 



The remaining terms here are the wave equation-like ones. We know that the solution 
is formed by two components (p -|- and p — </)), which propagate in opposite senses at 
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speed C (in general we need the matter part to be written in hyperbolic form, so it can be 
diagonalized) . 

Now here is where we need to put a boundary condition. What we want is: first, the outgoing 
component not to be touched; any condition over it would be superfluous or contradictory. 
Simply let it go out. Second, as we suppose that there is no scalar field outside the net, 
nothing is coming in (we disregard nonlinear effects between outgoing radiation and external 
metric), so the boundary condition to impose is simply 

Pi + ^1 = P2 + ^2 = (28) 

The third step is again S{^). 

This is all. The only boundary condition we need to put is (|28|), which is very simple 
indeed. Remember that there is not gravitational radiation in spherical symmetry. The 
momentum constraint plus the coordinate gauge condition fully determine the remaining 
(metrical) part of the boundary conditions. Note that the relation of D^g^ with holds 
when the momentum constraint is satisfied, so we can use the computed values of Vr at the 
boundary to obtain the value of the incoming non-gauge field there. 

VII. CONCLUSIONS 

The stability and accuracy near the origin of spherically symmetric dynamical systems in 
Numerical Relativity is strongly improved by a new formalism which provides an automatic 
care-free treatment of the geometrical factors arising from spherical symmetry. These treat- 
ment involves making use of the redundant variable previously introduced for other pur- 
poses (to ensure the hyperbolicity of the system). The same formalism also allows boundary 
treatment based on physical considerations and helps in increasing accuracy and stability. 
That is, the same characteristic of our system, i.e. the introduction of the new quantity 
Vri allows us to simultaneously improve the treatment of the two main problems in spher- 
ically symmetric star-like systems: the center (through a redefinition of which brings 
the system to a regular form) and the boundary (through hyperbolicity achieved thanks 
to which lets us distinguish between incoming and outcoming fields and treating them 
according to physical criteria). 

We have started from equations (ij-IH), which are just the ID version of a general 3D 
system which is actually being applied in multidimensional Numerical Relativity. Our 
point is that this system, when rearranged to deal with the geometrical factors, works fine 
in the ID case. 

We are presently considering the axially symmetric case (2D), where again the 'redundant' 
quantities Vi could be helpful to get rid of the geometrical factors that arise from the coor- 
dinate singularity at the symmetry axis. But this is still an open question. 
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FIGURES 



FIG. 1. The central value of the scalar field is graphed against time for an U-Branch star 
which has suffered a small perturbation (1% reduction on its central field). The star migrates from 
its former U-Branch position (central field value equal to 0.4) towards a S-Branch configuration 
corresponding to a central field value equal to 0.1, approximately. The migration consists in wide, 
diminishing oscillations of both the central field and the star radius around their final equilibrium 
configuration. 
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